% MT fwd & inverse computaton - keller & kaufman 
% latest date 04.01.2002
% inverse problems ----(???)
% fwd nlayer bug fixed 02.1.2002

% t = period (input - HZ)
% rho = layer resisitivity ohm-m
% h = layer thickness (IN KM)
% j9XY;
% Added junk for MCA students during  May-Juni 2005

mu = 4*pi*10^-7;
d1 = sqrt(1/mu);
%h   = [d1];
%rho = [1,0.01];
fid = fopen('jann_test.txt','wt');
frq = [0.0001 0.0003 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10 30 100 300 1000 3000 10000];
fprintf(fid,'%6s %6s %6s %6s %6s ','h1','h2','rho1','rho2','rho3');
    for j = 1:17,
        if frq(j) < 1,
            fprintf(fid,'%6.4f ',frq(j));
        elseif frq(j) >=1,
            fprintf(fid,'%6.0f ',frq(j));
        end;
    end;
  
    for j = 1:17,
        if frq(j) < 1,
            fprintf(fid,'%6.4f ',frq(j));
        elseif frq(j) >=1,
            fprintf(fid,'%6.0f ',frq(j));
        end;
    end;
    fprintf(fid,'tag - target');

fprintf(fid,'\n');
rhmin = log10(1);
rhmax = log10(6000);
phmin = 0;
phmax = 90;

t = 1./frq;
pk = 1;
for i = 1:1000,
    
    h1 = floor(rand*1000*2+100); %basalt layer
    h2 = floor(rand*1000*2+100); %sediment layer thickness
    rho1 = floor(rand*100+50);
    rho2 = floor(rand*20+5);
    rho3 = floor(rand*5000+1000);
    if rand <= 0.5,
        rho2=rho3; %no sediment;
        flg = 0;
    else,
        flg = 1;
    end;
    h = [h1 h2];
    rho = [rho1 rho2 rho3];
    [rh,ph]=mtfw(rho,h,t);
    
    rh = normrnd(rh,5); % gaussian error with sigma == 5
    ph = normrnd(ph,1);
    rh = rh/(31.5827340835/20); % Strange but this is a must
    
    rh = log10(rh);  % logarithm
    rh = (rh-rhmin)./(rhmax-rhmin); %normalizing between 0 and 1
    ph = (ph-phmin)./(phmax-phmin);%normalizing between 0 and 1


    %fprintf(fid,'%6d %6d %6d %6d %6d ',h1, h2,rho1 ,rho2 ,rho3);
    
    for j = 1:17,
        fprintf(fid,'%6.4f ',rh(j)); %normalized rho
    end;
    for j = 1:17,
        fprintf(fid,'%6.4f ',ph(j)); %normalized phase
    end;
       fprintf(fid,'%6d',flg); % flag for status
    fprintf(fid,'\n');
end;


fclose all;
